A nonlinear random walk approach to concentration-dependent contaminant 

transport in porous media 
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We propose a nonlinear random walk model to describe the dynamics of dense contaminant 
plumes in porous media. A coupling between concentration and velocity fields is found, so that 
transport displays non-Fickian features. The qualitative behavior of the pollutant spatial profiles 
and moments is explored with the help of Monte Carlo simulation, within a Continuous Time 
Random Walk approach. Model outcomes are then compared with experimental measurements of 
variable-density contaminant transport in homogeneous and saturated vertical columns. 



I. INTRODUCTION 

Non-Fickian (anomalous) transport is a widespread 
feature of contaminant migration in porous media 
Specifically, 'non-Fickian' means that the spread of the 
transported species grows nonlinearly in time, {x^it) — 
{x{t))'^) ^ t^, (3 ^ 1, the resulting concentration profiles 
displaying a non-Gaussian behavior [J, |3i 0] ■ This is in 
contrast with the linear spread and Gaussian shapes usu- 
ally expected for particles migration in perfectly homo- 
geneous media, where the Fickian advection-dispersion 
equation apphes: see, e.g., Q and References therein. A 
broad spectrum of physical reasons have been invoked 
to explain the observed deviations from Gaussianity. 
For instance, the homogeneity hypothesis becomes ques- 
tionable in presence of irregularities at multiple space 
scales d, 0, complex structures of flow streams 0, Q 
and saturation distribution within the medium and 
physico-chemical exchanges of the pollutant particles 
with the surrounding material Another important 

source of non-Fickian behaviors is the collective motion of 
pollutants due to reciprocal interactions. A well-known 
example is provided by reactive transport, where two or 
more chemical species may combine (reversibly or irre- 
versibly) to give birth to new ones. Even in homoge- 
neous media, this may lead to intricate contaminant pat- 
terns Rd|, whose complexity could be further increased 
by the presence of spatial heterogeneities [12]. 

Intuitively, the dynamics of concentrated particles will 
also display nonlinear, collective phenomena. Indeed, 
the motion of a single pollutant parcel depends on the 
density of the surrounding fluid, which in turn is af- 
fected by the number of such parcels nearby, so that 
the microscopic trajectories are correlated. Transport 
of dense pollutant plumes has been long investigated, 
yet kee^s raising many conceptual as well a,s practical 
issues [il, lii. Ill III O E H H ill, [ll. stud- 
ies cover both homogeneous saturated and heterogenous 
unsaturated materials [2^ [25l [26l . [13] : extensive reviews 
may be found, e.g., in j28l. |29|T^ Strong density gradi- 
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ents are encountered when either the contaminant itself 
is highly concentrated at the source, or the plume flows 
through regions that are rich in salt; in particular, this 
latter case might become a major concern for radioactive 
waste disposal near salt domes JiOj ■ 

Similarly as Brownian motion is related to the diffu- 
sion equation, concentration-dependent particles paths 
can be formally shown to lead to a family of nonlin- 
ear Fokker-Planck transport equations, on the grounds 
of a statistical-mechanical approach; see, e.g., |3ll [3^ 
for a detailed account of recent advances. The displace- 
ments of a particle in the medium are thought to be 
affected by the number of other particles in its initial 
or final position, or both [33] : this allows better un- 
derstanding the small-scale dynamics, rather than im- 
posin g th e macrosco pic eg nations on a phenomenological 

basis [Hi [si, [si, [33, lal. 

Adopting a somewhat similar perspective, we propose 
here a simple model for the collective concentration- 
dependent dynamics of a dense contaminant plume and 
explore its qualitative behavior by resorting to Monte 
Carlo simulation. Model predictions are then validated 
on experimental data. This paper is organized as follows: 
in Section ini we develop a stochastic equation that de- 
scribes the motion of a pollutant parcel in a dense fluid. 
In Section IIIIl we discuss the qualitative behavior of the 
model and the interplay of its components. Then, in Sec- 
tion llVl we proceed to compare model outcomes to exper- 
imental results of variable-density contaminant transport 
in saturated homogeneous porous columns. Finally, the 
potentialities and the limits of the proposed approach are 
evidenced in Section (Vl 



II. A NONLINEAR TRANSPORT MODEL 

Let us consider a vertical column filled with sand. For 
sake of simplicity, we start by assuming that the sand 
is uniformly packed and well-mixed, so that the porous 
medium can be considered as homogeneous, and that the 
column is fully saturated in water. When the ratio be- 
tween the length and the diameter of the column (the so- 
called aspect ratio) is much greater than one, the system 
can be regarded as one-dimensional, to a first approxima- 
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FIG. 1: Concentration profiles at time t = 0.45 h, for step 
injection from t = h to t — 0.23 h. Solid line represents 
Fickian transport (e = 0); dotted line (injection from the 
top) and dashed line (injection from the bottom) represent 
nonlinear concentration-dependent transport. 



FIG. 2: Variance of the particles plume as a function of time, 
for step injection from t = Ohtot = 0.23 h. Fickian (dots, e = 
0) and concentration-dependent (dotted line: injection from 
the top; dashed line: injection from the bottom) transport 
processes are displayed. 



tion. Suppose now that a given amount of contaminant 
fluid is injected into the column: we can conceptually rep- 
resent the pollutant plume as a collection of fluid parcels 
i — 1, N, each containing a fraction rrii = M/N of the 
total contaminant mass M . When the effects of molecu- 
lar diffusion are negligible, it is reasonable to assume that 
rrii will not change in the course of plume evolution [s^ . 
If V is the reference volume of the injected pollutant, 
each parcel carries a volume Vi = V/N. 

The projections of forces acting on a parcel i in the di- 
rection of the flow are: the pressure gradient imposed by 
the injecting pump, Fp, supposedly constant; the viscous 
resistence which opposes flow, namely Fy = —jUi(t), 
where the friction coefficient 7 = /i/fc is given by the 
ratio of the fluid dynamic viscosity /i [Kg/m s] and the 
medium permeability k [m^], and Ui(t) is the local ve- 
locity of a parcel; gravity and buoyancy, which can be 
written as Fg = g{pi — pj), where g is the gravity ac- 
celeration. Pi is the density of the contaminant parcel 
and p( is the density of the fluid surrounding the par- 
cel i. Mechanical dispersion can be taken into account 
by adding stochastic fluctuations Si around the parcel 
velocity Ui{t) ,37]. It is customary to assume 



Si cx ^y\{ui)\rii 



(1) 



where (ui) is the ensemble average of the particles veloc- 
ities (provided that the medium is sufflciently homoge- 
neous [38] ) and i]i is an uncorrelated white noise with zero 
mean and unit variance. The constant of proportionality 
determines the strength of the velocity fluctuations and 
is thus related to the dispersivity a [m] of the porous 
material. Then, the forces balance reads 



PiUi 



F„ 



-fu, + g{p^~ pj) + S^, 



(2) 



where the reference axes system is chosen so that gravity 
is positive pointing downwards and the explicit depen- 
dence on time has been omitted. It appears that the 
absolute value of the pollutant density does not play a 
major role, the plume migration being mostly controlled 
by relative density differences: this is coherent with ex- 
perimental evidences [H Hi, [13, El, IH, El . 

Let us now focus on the two terms pi and p{ . The 
density of a contaminant parcel can be expressed as 



Pi = 



PoVi 



(3) 



where poVi is the mass of reference fluid (e.g., water) 
contained in Vi and po its density. By resorting to the 
deflnitions of to,- and w, , we obtain 



1 M 
Po V 



(4) 



where finally M/V is given by the product of the mo- 
lar concentration (7™°' [mol/L] times the molar mass 
[g/mol] of the injected species. Even modest density dif- 
ferences with respect to the resident fiuid (of the order 
of a few percents) might sensibly affect the contaminant 
dynamics [H, [33|. Hence, we focus on this case and 
think of Pi as a small perturbation compared to po, i.e., 
e = {M/V)/po < 1. As for the local fluid density p{ , 



J - 



PQdx + m{xi, t) 
dx ' 



(5) 



where m{xi,t) is the pollutant mass contained in an ele- 
mentary volume dx around the position Xi of the parcel i. 
We are assuming that each parcel is aware of the presence 
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FIG. 3: Concentration profiles at time t = 0.48 h, for step 
injection from t = h to t = 0.24 h. Solid line represents 
anomalous transport due to spatial heterogeneities, modelled 
by a waiting times pdf with power-law decay iP{t) ~ f-'^l"^ 
(with £ = 0); dotted line (injection from the top) and 
dashed line (injection from the bottom) represent nonlinear 
concentration-dependent transport coupled with the effects of 
the spatial heterogeneities, for the same iP{t). 
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FIG. 4: Variance of the particles plume as a function of time, 
for step injection from t = Ohtot = 0.24 h. Anomalous trans- 
port due to spatial heterogeneities, modelled by a waiting 
times pdf with power-law decay i/'(r) ~ t-~3/2 ^^j^Jj ^ _ g-j jg 
represented with a solid line. Concentration-dependent trans- 
port, coupled with the effects of the spatial heterogeneities, is 
displayed as dotted line (injection from the top) and dashed 
line (injection from the bottom), for the same i/'(r). 



of the others only at short range, through the effects of lo- 
cal density variations. Since m{xi,t) = n{xi,t)mi, where 
n{xi,t) is the number of pollutant parcels in [xi,Xi + dx] 
at time t, we can finally rewrite 



PI 



Po 



1 



n{xi,t) 

Nn 



(6) 



where A'o is a dimensionless normalization factor such 
that M/V — Norrii/dx. In practice, iVo expresses the 
(arbitrary, but sufficiently large) number of contaminant 
parcels that are initially attributed to each dx to rep- 
resent the average density AI/V at injection. At each 
time step, the quantity c{xi,t) = n{xi,t)/NQ identifies 
the contaminant concentration at position Xi. 

The role of viscosity has been condensed in the con- 
stant parameter 7. In reality, viscosity depends on con- 
taminant concentration, but its variations are frequently 
less relevant than those of density and are thus ne- 
glected [ll]. Within the proposed formulation, in- 
cluding a functional dependence of the kind 7 = 7o(l + 
Xc{xi, t)), where 70 is the reference value in the fluid and 
A is a (small) constant, would be straightforward. In 
the following, however, we always suppose that 7 ~ 79. 
Moreover, we do not address the possible dependence of 
density and viscosity on other physical variables, such as 
temperature. 

Finally, assuming that inertial effects can be neglected 
(which is the case, provided that viscous forces are dom- 
inant), and making use of expressions [5] and [51 we can 
rewrite Eq. [2] in Langevin form 



Equation [7] describes the random walk of a fluid parcel 
which is advected at a concentration-dependent speed 
u{c) ^ Up + Ug + Uc, with Up = j~^Fp, Ug = j^'^gpoe 
and Uc — —j^^gpQec{xi,t), and dispersed by fluctuations 
whose amplitude is std{"f~^ Sidt) = [2a\{ui{t))\dt]^^^ . 
Note that dispersion D{c) — a\{ui{t))\ is also a func- 
tion of concentration, through the dependence on the 
ensemble-averaged velocity. 

By relying upon the results resumed in, e.g., (33 |. it 
is possible to show that the smoothed contaminant con- 
centration held c{x,t) — 5{x — Xi{t))) corresponding 
to particles undergoing the random walk in [7] obeys a 
nonlinear Fokkcr-Planck equation 



d-&^'^ = -d-x 



u(c(x,t)) - ^D{c{x,t)) 



c(x,t). (8) 



Xi = u{c) -f 7 ^Si 



(7) 



Equations of this form are well-known and com- 
monly arise in the context of transport processes with 
concentration-dependent dispersion and/or velocity: see, 
e.g., While we will not make explicit use of its prop- 
erties in the following, Eq. [5] provides the necessary link 
between the microscopic stochastic particles dynamics in 
Eq. [7] and the deterministic evolution of the associated 
ensemble-averaged concentration field. Note that the ef- 
fects of mutual interactions in Eq.[H]bccomc negligible for 
e — > 0, i.e., when the molar concentration of the injected 
solution is weak. In this case, the particles trajectories 
are independent, Uiit) Up, and Eq. [8] degenerates to a 
standard advection-dispersion equation, so that Fickian 
transport is recovered, with D = aup. For a given value 
of e > 0, the nonlinear coupling plays a minor role at 
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FIG. 5: Downwards injection at a reference molarity (7™°' = 
0.2 mol/L. Contaminant concentration curves ce{t) measured 
at sections £ = 7.7,23.1,38.5,46.2, and 77 cm (from left to 
right), as a function of time. Squares correspond to experi- 
mental data, solid lines to Monte Carlo simulation. 

short time scales also for n{xi, t) —> 0, i.e., when the num- 
ber of contaminant particles in the considered dx is small. 
This is the case when dispersion dominates, so that fluid 
parcels are rapidly dragged far apart and can hardly in- 
teract. Eventually, at longer time scales, dispersion will 
usually overcome the effects due to concentration. 

III. DISCUSSION 

Equation [7] defines a discrete-time random walk where 
the particles positions are updated at each time step dt. 
In view of the possibility of describing a broad class of 
porous materials, such as heterogeneous and/or unsatu- 
rated media, it is expedient to resort to the more gen- 
eral Continuous Time Random Walk (CTRW) formal- 
ism [lll45|, where particles trajectories alternate random 
jumps (drawn from a pdf p(s)) and random waiting times 
(drawn from a pdf i^ir)) at each visited spatial site. The 
pdf Tpir) identifies the velocity spectrum in the traversed 
material: flows in homogeneous porous media such as 
those considered here (where it is reasonable to assume 
that the sojourn times at each site must be on average the 
same [l|) correspond to choosing a Poisson pdf tpir), so 
that a single time-scale, e.g., the average (r) of the distri- 
bution, dominates [ij. As for the displacements, the spa- 
tial scales of advection and dispersion in the CTRW are 
determined by the cumulants k of the jump lengths dis- 
tribution p(s) and are not separated a priori [Ij. A com- 
mon choice is to adopt a Gaussian pdf p(s), so that the 
first two cumulants are sufhcient to characterize trans- 
port: in particular, ni is associated to advection and K2 
to dispersion. 

We can then rephrase the stochastic process defined in 
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FIG. 6: Downwards injection at a reference molarity (7™° = 

0. 2 mol/L. Moments (t*(^)) of passage times t{£\ as a func- 
tion of various column heights I. Crosses represent the mean 
of the passage times (fc = 1), circles the second moment 
(fe = 2); the latter has been divided by a factor of 20 in 
order to have comparable scales. Solid (k = 1) and dashed 
lines (fc = 2) are the results of Monte Carlo simulation. 

Eq. [7] by resorting to a CTRW where waiting times obey 
a Poisson pdf with mean (r), i.e., 

^(r) = 7\e-/<-) (9) 

and jumps obey a Gaussian pdf with concentration- 
dependent cumulants k\ = u{c){t) and K2 = 2D{c){t), 

1. e., 

p{s) - p{s\c) = . (10) 

The CTRW formalism has been here introduced on phe- 
nomenological basis, as a generalization of Eq.[71 hints for 
a rigorous derivation of concentration-dependent transi- 
tion rates within a nonlinear master equation formula- 
tion are provided, e.g., in The process defined by [9] 
and 1101 can be easily simulated by Monte Carlo method. 
Note however that the jump lengths distribution p{s\c) 
explicitly depends on concentration, so that particles tra- 
jectories are not independent and mutual interactions re- 
quire knowing the concentration field (i.e., the locations 
of the entire ensemble) before updating walkers positions. 
Starting from a known initial condition c(a;,0), particles 
are displaced at each time step by drawing waiting times 
and jump lengths from [5] and 1101 respectively, and the 
new concentration field is recursively determined for the 
following time step. Walkers whose random waiting time 
is longer than the current time step stay in the same spa- 
tial site. For continuous contaminant spills, new particles 
are added at the column entrance for the duration of the 
injection. In order to attain convergence in simulations, 
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FIG. 7: Upwards injection at a reference molarity C"^"' —0.1 
mol/L. Contaminant concentration curves Ci{t) measured at 
sections £ = 7.7, 23.1, 38.5, 46.2, and 77 cm (from left to right), 
as a function of time. Squares correspond to experimental 
data, solid lines to Monte Carlo simulation. 



we could either i) choose very small time steps for up- 
dating displacements and concentration field, or ii) at 
each (larger) time step iteratevely compute the values of 
displacements and concentration until their relative error 
is below a given threshold. After testing both methods, 
we found convenient to resort to the former: the optimal 
value of the time step was determined by trial-and-error. 

Monte Carlo simulation offers an expedient means of 
exploring the qualitative features of the nonlinear CTRW 
transport model described above. In particular, we pro- 
ceed now to analyze spatial contaminant concentration 
profiles (at fixed time) for small values of the parameter e; 
this allow getting insights on the relevance of the coupling 
between velocity and concentration. Figure [T] compares 
a Fickian contaminant profile, corresponding to e = 0, 
with typical spatial profiles for concentration-dependent 
transport (0 < e <C 1). Downwards injection gives rise 
to positively skewed profiles, while the opposite is true 
for upwards injection. In all cases, we considered a step 
injection of finite duration. Nonlinear transport clearly 
displays asymmetric profiles, whereas Fickian transport 
corresponds to Gaussian (symmetric) profiles. 

The time duration of the contaminant injection is a 
key factor in determining the spatial shape of the plume. 
Indeed, the coupling between concentration and veloc- 
ity is in competition with dispersion, which in turn is 
induced by the average velocity {ui{t)). The stronger 
the velocity Up, the lesser is the relevance of the nonlin- 
ear term in u{c). The injected plume might have such 
a limited extension that dispersion rapidly dominates 
concentration-dependent effects: in other words, because 
of their velocity, fluid parcels become quickly dispersed, 
and their interactions through the density field are weak. 
This prediction is coherent with our experimental mea- 




20 40 60 80 



e \cm\ 

FIG. 8: Upwards injection at a reference molarity (7™°' =0.1 
mol/L. Moments {t^{(-)) of passage times t{t), as a function 
of various column heights I. Crosses represent the mean of 
the passage times [k = 1), circles the second moment {k = 2); 
the latter has been divided by a factor of 20 in order to have 
comparable scales. Solid {k — 1) and dashed lines {k = 2) are 
the results of Monte Carlo simulation. 



sures: increasing the imposed flux (at flxed C™°'), the 
contaminant proflles approach standard Fickian shapes. 
At the opposite, the longer the extension of the injected 
plume and the more persistent are the effects of the re- 
ciprocal interactions, before eventually dispersion takes 
over. This phenomenon has already been experimentally 
detected for the case of viscosity-dependent transport of 
'slices' of finite duration [46] . Therefore, for a given value 
of e, the relevance of the nonlinear coupling is stronger 
for small velocity fields Up and long injection times. 

On the basis of the observations above, one might ex- 
pect that the effects of the nonlinear coupling would come 
into play mainly through velocity variations. Actually, it 
turns out that the average particle velocity {ui{t)) is only 
slightly affected by density, provided that e is not too 
large. On the other hand, fluctuations around the mean 
velocity (induced by the nonlinear terms) do not sim- 
ply average out, and contribute instead to an apparent 
plume dispersion (in addition to a\{ui{t))\ ~ aup). This 
is a relevant and subtle outcome, which is ultimately re- 
sponsible for the skewed shape of the pollutant proflles, 
the concentration-dependent contribution to dispersion 
being proportional to density differences (and thus non- 
symmetric) . 

In Fig. [2] we display the behavior of the contaminant 
variance (x^(t) — {x{t))'^) as a function of time, as com- 
puted by Monte Carlo simulation. For the case of Fick- 
ian transport (e=0), the variance is a straight line, as 
expected. For concentration-dependent transport, the 
variance turns out to be a nonlinear function of time and 
appreciably deviates from the Fickian behavior. This is 
indeed the hallmark of anomalous diffusion. On the con- 
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IV. COMPARISON WITH EXPERIMENTAL 
RESULTS 
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FIG. 9: Skewness Xiil) and kurtosis X4(^) of the arrival times 
considered in Figs. [5] and [T] as a function of column height I. 
Estimates of X'a{^) for experimental data: squares (upwards 
injection) and crosses (downwards); Monte Carlo estimates: 
dashed lines. Estimates of x^i^) for experimental data: circles 
(upwards injection) and triangles (downwards); Monte Carlo 
estimates: solid lines. Estimates for Gaussian transport (e = 
0): dotted {Xi{^)) ^-nd dotted-dashed {Xi{^)) lines. 



trary, the average of the contaminant plume (not shown 
here) is found to be linear in time, for small values of e. 

While in the present discussion we have made the as- 
sumption of considering flows in homogeneous porous 
media, which amounts to drawing waiting times from a 
Poisson pdf, the CTRW framework straightforwardly al- 
lows taking into account spatial heterogeneities, due, e.g., 
to different grain sizes or variable saturation. The broad 
velocities spectra that are commonly found in hetero- 
geneous and/or unsaturated materials are mirrored in a 
broad distribution of time scales for the jumping rates be- 
tween sites: it is customary to incorporate such physical 
processes in iP{t) by considering power-law waiting times 
between consecutive displacements, possibly with an ex- 
ponential cut-off [l], 0, [ll]. Particles trajectories would 
then be affected on one hand by concentration-dependent 
displacements and on the other hand by anomalously 
long sojourns: the two processes may superpose, de- 
pending on the respective time scales. The qualitative 
behavior of the competition between density effects and 
heterogeneities is displayed in Figs. [3] and [H where we 
show spatial contaminant profiles and particles variance, 
respectively, for a waiting times pdf ')/'(''') ~ t""^/^. In 
particular, we remark that the asymmetry that was ev- 
ident for homogeneous transport (Fig. [J) is now hidden 
by the long tails of the pollutant profiles. 



In this Section, we test the proposed random walk 
model on some experimental measurements of dense con- 
taminant transport obtained at the Physical-Chemistry 
Department (DPC), CEA/Saclay. The experimental de- 
vice, named BEETI, consists of a dichromatic X-ray 
source (20 — 40 keV, 50 — 75 keV), applied to a verti- 
cal column of height H = SO cm and diameter D — 5 
cm (the aspect ratio is therefore H/D = 16 ^ 1). The 
X-ray transmitted countings allow quantitatively assess- 
ing the contaminant concentration inside the column (as 
a function of time), at various sections i: we denote this 
quantity by Ci{t). The different positions are explored by 
means of a remotely controlled rack rail that displaces the 
X-ray emitter and the coupled Nal detector. At the exit 
of the column, cg^nit) coincides with the breakthrough 
curve, which is the most frequently measured variable in 
contaminant migration experiments In the specific 
context of dense contaminant transport, only a few works 
have investigated the behavior of breakthrough curves 
corresponding to finite-duration injections, whereas at- 
tention is usually focused on the mixing properties at 
the interface between two layers of semi-infinite exten- 
sion (see, e.g., [iBIlll and References therein). 

The BEETI experimental setup allows for downwards 
as well as upwards fluid injection, and several kinds of 
flow regimes and porous materials can be tested, at var- 
ious saturation and/or heterogeneity conditions. To set 
the ideas, in the following we refer to fully saturated 
columns filled with homogeneously mixed Fontainebleau 
sand (bulk density 1.77 ± 0.01 g/cm^), with average 
grain diameter 200 /im. The average porosity is 6* = 
0.333 ±0.005 and the dispersivity is a = 0.1 cm. The ref- 
erence saturating fluid is water containing dissolved KCl 
(molar mass equal to 74.5 g/mol) at a molar concentra- 
tion of 10-3 mol/L, so that po = 998.3 Kg/m^ at T = 20 
C°. The injected contaminant is KI (molar mass equal to 
166 g/mol), at different molar concentrations. All mea- 
surements are performed at constant room temperature 
T = 20 C°. We estimated - 5 cm/h. Con- 

taminant flow is imposed at one end of the column and 
collected at the other end, where an electric conductivity 
meter provides a supplementary (independent) measure- 
ment of the breakthtrough curve. The pump imposes a 
steady state Darcy flow of g = Upd = 2 cm/h, which is 
verified by weighing the outgoing solution. The experi- 
mental conditions are such that clogging or formation of 
colloidal particles, which could alter the interpretation 
of the obtained results, can be excluded. Chemical re- 
actions or sorption/desorption phenomena can be ruled 
out as well. 

A representative example is shown in Figs. [5] and [5] 
for downwards injection of KI at g = 2 cm/h, with 
(jrnoi ^ 0.2 mol/L, so that e = 0.033. The time du- 
ration of injection is 3 h. Figure [H] compares the exper- 
imental concentration profiles (squares) with the Monte 
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Carlo simulation results (solid lines). From the point of 
view of Monte Carlo simulation, the quantity ci(t) is es- 
timated as the number of particles that are contained 
in a volume dx around the position £, at a given time 
t. In other words, ci(t) represents the distribution of 
the passage times at fixed positions. In principle, knowl- 
edge of the physical constants completely determines the 
free parameters of the simulation; in practice, however, a 
trial-and-error fine fitting around e and a is required in 
order to account for uncertainties. Despite the many as- 
sumptions and simplifications introduced in the random 
walk model, a good agreement is found between simula- 
tion and data. This agreement, moreover, is preserved 
all along the measurement points £, thus meaning that 
the proposed model allows capturing the full spatial dy- 
namics of the plume. It is evident that the asymmetric 
spatial shape that had been predicted on the basis of 
random walk simulations (Fig. [1]) is now mirrored in the 
shape of ci{t). Due to the interplay of concentration and 
velocity, a part of the contaminant plume is descending 
faster than the bulk. 

The agreement between model and experimental data 
is further substantiated by Fig. [6l where we compare 
the first two moments (t'^(£)), k = 1,2, of the passage 
times t{£) along the column. We remark that the slope 
of {t^{tj) is very close to the value which is con- 

sistent with the average particles velocity being almost 
unaffected by the concentration field. These findings 
are coherent with experimental observations and mod- 
els of density-dependent transport proposed in litera- 
ture [3, m, SO, SH El El. 

Comparable results have been obtained also for up- 
wards injection. A representative example is shown in 
Figs. [7] and [H for q^2 cm/h and C""°' = 0.1 mol/L, so 
that e = 0.017. The time duration of injection is 3 h. 
The asymmetric tail of the contaminant concentration 
profiles is now on the right, meaning that part of the 
bulk is delayed because of density effects (cf. Fig. [7]). A 
slightly less satisfactory agreement is found for the pro- 
files at intermediate heights, which could be attributed 
to neglecting inertial contributions in Eq. [2] Nonethe- 
less, the breakthrough curve and the moments (Fig. [5]) 
are well captured by the random walk model. 

Finally, in order to emphasize the departure of the con- 
centration profiles shown in Fig. [5] and [7] from Gaussian 
behavior, in Fig. [5] we provide the skewness Xi{^) ^-i^d 
kurtosis Xi{^) of t^i^ arrival times El, as a function of 
column height I. Experimental data estimates lie close 
to those of Monte Carlo simulations. For comparison, 
the case of Gaussian transport (i.e., e = 0) is plotted 
in the same figure: the difference with respect to dense 
contaminant transport is clearly noticeable. Remark in 
particular that for downwards injection Xsi^) < and de- 
creases with £, whereas for upwards injection Xsi^) > 
and increases with £. For Xi{^)j similar deviations from 
Gaussian behavior are observed, though partially hidden 
by limited statistics. 

In principle, one might wonder whether a standard lin- 



ear CTRW with algebraic ■!/'(^), which also gives rise to 
asymmetric breakthrough curves with long tails, could be 
applied to fit the experimental data. However, this hy- 
pothesis is in contrast with two basic facts: first, adopt- 
ing a power-law waiting time pdf is somehow unjustified, 
since the medium is homogeneous; second, a standard 
linear CTRW approach could not explain why the asym- 
metry of the breakthrough curves is affected by the flow 
direction. So far, our experimental activities have ex- 
clusively concerned the transport of dense contaminant 
plumes in homogeneous saturated columns. However, 
further tests are in order, to explore the case of hetero- 
geneous and/or unsaturated porous media. The BEETI 
device, thanks to the dual-energy source, can determine 
at the same time contaminant concentration and water 
content at each section: it would be thus interesting to 
compare model predictions (Figs. [3] and |4]) with experi- 
mental data. 



V. CONCLUSIONS 

We have proposed a nonlinear random walk approach 
to the modelling of variable-density contaminant flows in 
porous media, within a CTRW framework. The qualita- 
tive behavior of this model has been explored by means 
of Monte Carlo simulation: particles trajectories are cor- 
related via the density fleld, so that transport is non- 
Fickian and the plume variance grows nonlinearly in 
time. When the molar concentration of the injected pol- 
lutants is similar to that of the resident fluid, the usual 
Fickian behavior is recovered. Within CTRW, it is pos- 
sible to describe transport through both homogeneous 
and heterogeneous materials: in this latter case, we have 
shown that the effects of concentration-dependent dy- 
namics are in competition with (and might partially hid- 
den by) those of spatial heterogeneities. 

The proposed random walk model is admittedly sim- 
ple, since the full spectrum of interactions that actually 
take place between the velocity and density fields has 
been condensed in a single nonlinear coupling at the scale 
of particles trajectories. Detailed studies show that the 
physics behind variable-density transport is essentially 
3d, or at least 2d, because of the complex interfacial dy- 
namics between two fluids of different densities and/or 
viscosities [H [H 113, El, [11 HI HI, HI . Neglecting these 
phenomena leads to descriptions that must be necessar- 
ily intended in a mean-field sense: only the coarse-scale 
behavior of the real system can be captured, and the fine- 
scale details are averaged out [H io. El, Ell E3, El, E3 ■ 
Moreover, we have made the hypothesis that molecular 
diffusion is negligible with respect to mechanical disper- 
sion, and that viscosity can be considered as constant, to 
a first approximation. 

Yet, our random walk model compares well to a set 
of dense contaminant transport measurements realized 
by means of the BEETI device. The experimental con- 
ditions ensure that most of the introduced simplifica- 
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tions actually apply: the aspect ratio of the column is 
large, so that migration is almost Id; viscosity varia- 
tions are weaker than density variations; molecular dif- 
fusion is smaller than dispersion. It seems reasonable to 
think that the limits of validity of the proposed model 
will clearly emerge when these hypotheses are not veri- 
fied: experimental activities are ongoing and will be pre- 
sented in a forthcoming work. In particular, we expect 
our model to provide a satisfactory agreement with mea- 
sured data when it is possible to consider density varia- 
tions as small perturbations with respect to the resident 
fluid (i.e., e <C 1). For larger density differences, other, 
more complex couplings should perhaps be introduced, 
possibly involving higher-order nonlinearities and long- 
range correlations. The findings in [so], for instance, sug- 
gest that in presence of relevant density gradients even 
the validity of Fick and Darcy laws at microscopic scale 
should be carefully reconsidered. In this respect, Monte 
Carlo simulation might be complemented, e.g., by the 
promising computational tool of Smoothed Particle Hy- 



drodynamics, which has been recently applied with suc- 
cess to the numerical study of variable-density flows with 
stochastic dispersion [23 |. 

The proposed random walk approach has been moti- 
vated by a specific problem in contaminant migration; 
many other physical processes where the CTRW for- 
malism applies may exhibit particles paths correlated 
via the density field, so that relevant advances would 
be achieved by formally generalizing the CTRW the- 
ory for the case of concentration-dependent distributions. 
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